A previous presentation considered numerical visualization of the eigenvalues and eigenfunctions of a symmetric one-dimensional power potential with arbitrary real exponent. The Schrödinger equation that was employed is

$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}+\frac{m{\omega }_a^2}{2}|x{|}^a\psi =E\psi$

where the potential has a (nonanalytic) absolute value for symmetry about the origin, and the subscript on the frequency is an indication that its dimensioned units must change to accommodate the difference in spatial behavior from the simple harmonic oscillator. The visualization was performed by expanding the eigenfunctions for the potential with arbitrary exponent in terms of eigenfunctions for the quantum mechanical simple harmonic oscillator, since the latter are relatively simple mathematically and readily amenable to numerical evaluation.

As pointed out in that presentation, the simple harmonic oscillator is one of a handful of potentials with full analytic solutions to the Schrödinger equation, and along with the mathematically related Coulomb potential has one of the simplest solutions. The previous numerical evaluations are useful for having an idea of the behavior of eigenvalues and eigenfunctions during alteration of the exponent of the power potential, but what one really wants is something closer to the analytic solution for the simple harmonic oscillator. This can be achieved with an application of WKB theory to the bound states of the power potential with arbitrary exponent.

This presentation will make use of methods developed in Chapter 10 of Bender & Orszag, Advanced Mathematical Methods for Scientists and Engineers, a delightful compendium of many a useful technique. It will give an approximate analytic form to the numerical visualizations of the values and functions found in the previous presentation, and extend the results from power potentials with positive exponents to those with exponents greater than minus two.

The Schrödinger equation as given has a singularity at infinity with an s-rank of $\frac{a+4}{2}$ (see this discussion). As long as $a>-2$ this singularity is irregular, and the standard way of approximating asymptotic behavior at an irregular singularity is with an exponential. The general method of WKB theory is to expand the function in the exponential in powers of Planck’s constant:

$\psi \left(x\right)=exp\left[\frac{1}{\hslash }\sum _{k=0}^{\infty }{\hslash }^k{S}_k\right]$

The series in the exponent is in general asymptotic: the first few terms will give a decent approximation to the analytic solution, but inclusion of too many terms leads to a function that diverges from the exact solution. The second term S₁ is a logarithm as will be seen shortly.

To apply the method, first write a general one-dimensional Schrödinger equation in normal form,

${\hslash }^2\frac{d^2\psi }{d{x}^2}=Q\left(x\right)\psi Q\left(x\right)=2m(V-E)$

then insert the expansion and evaluate the second derivative:

$\sum _{k=0}^{\infty }{\hslash }^{k+1}{S}_k^{″}+[\sum _{k=0}^{\infty }{\hslash }^k{S}_k^{\prime }{]}^2=Q$

Collecting on powers of $\hslash$, the terms in the expansion are determined by

$\begin{array}{c}(S_0^{\prime }{)}^2=Q\\ \\ S_{k-1}^{″}+\sum _{l=0}^k{S}_l^{\prime }{S}_{k-l}^{\prime }=0,k\ge 1\end{array}$

This part of the WKB process is essentially an application of the method of dominant balance, where the inclusion of Planck’s constant makes the determination straightforward. The second term of the series is determined solely by the first,

$\begin{array}{c}{S}_0^{″}+2S_1^{\prime }{S}_0^{\prime }=0\\ \\ S_1=ln(S_0^{\prime }{)}^{-1/2}=ln{Q}^{-1/4}\end{array}$

which is the logarithm mentioned. For purposes of this presentation, only the first two terms in the series will be retained, so that the approximation to the exact wave function is proportional to either one of the two functions

$\psi \left(x\right)\approx Q^{-1/4}exp[\pm \frac{1}{\hslash }\int dx\sqrt{Q}]$

where the limits of integration must be specified explicitly in order to determine eigenvalues.

Inside a potential well energy is greater than the value of the potential, so that Q < 0 and the argument of the exponential is imaginary. This means that the wave function can be represented by circular functions in the potential well. Outside of the well energy Q > 0 and one can choose the decreasing exponentials for a physically meaningful bounded wave function.

The approximation has a problem, however, at turning points of the classical potential where V = E because the first inverse factor becomes infinite at those points. The breakdown in the approximation occurs because the differential equation changes radically at turning points: the right-hand side of the equation disappears. To estimate the behavior through turning points, replace the function Q by a linear term with a constant coefficient, which is the first power that appears in a power-series expansion of the function. The differential equation near turning points is then

${\hslash }^2\frac{d^2{\psi }_{\mathrm{turn}}}{d{x}^2}\approx c_{\mathrm{turn}}(x-x_{\mathrm{turn}}){\psi }_{\mathrm{turn}}$

This is the Airy differential equation, and that means that the wave function on either side of the turning point can be estimated with asymptotic forms of Airy functions. Justification for these estimations are given fully in Bender & Orszag and will not be considered here. Wave function approximations on either side of turning points are then related by so-called connection formulas between asymptotic forms of Airy functions.

In the region to the right of the potential well, the wave function will be taken to be

${\psi }_{\mathrm{right}}\left(x\right)=\frac{1}{2}{c}_{\mathrm{right}}{Q}^{-1/4}exp[-\frac{1}{\hslash }\underset{x_{\mathrm{turn}}}{\overset{x}{\int }}dx\sqrt{Q}]$

Moving inside the potential well, the connection formula for an Airy function of the first kind alters the exponential part of this wave function to

${\psi }_{\mathrm{inside}}\left(x\right)=c_{\mathrm{right}}(-Q{)}^{-1/4}sin[\frac{1}{\hslash }\underset{x}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}+\frac{\pi }{4}]$

The sign under the radial accounts for the explicit sign change of Q when moving into the potential, which is why the exponential becomes a circular function. The reversal of the order of the limits of integration represents the negative argument for the Airy function in the region to the left of the turning point, and the fractional phase shift comes from the Stokes constant appropriate to this connection.

In the region around the left turning point, the linear approximation to the integrand changes sign, and this needs to be reflected in the approximation to the wave function. The wave function to the left of the potential well is thus

${\psi }_{\mathrm{left}}\left(x\right)=\frac{1}{2}{c}_{\mathrm{left}}{Q}^{-1/4}exp[-\frac{1}{\hslash }\underset{x}{\overset{-x_{\mathrm{turn}}}{\int }}dx\sqrt{Q}]$

Moving inside the potential well, the connection formula for an Airy function of the first kind alters the exponential part of this wave function to

${\psi }_{\mathrm{inside}}\left(x\right)=c_{\mathrm{left}}(-Q{)}^{-1/4}sin[\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x}{\int }}dx\sqrt{-Q}+\frac{\pi }{4}]$

And now for some fun! The quantization condition for allowed energy values arises from ensuring the the two forms for the wave function inside the potential are identical. Rewriting this second form,

$\begin{array}{l}{\psi }_{\mathrm{inside}}\left(x\right)=c_{\mathrm{left}}(-Q{)}^{-1/4}sin[\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}-\frac{1}{\hslash }\underset{x}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}+\frac{\pi }{4}]\\ \phantom{{\psi }_{\mathrm{inside}}\left(x\right)}=-c_{\mathrm{left}}(-Q{)}^{-1/4}sin[\frac{1}{\hslash }\underset{x}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}-\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}-\frac{\pi }{4}]\\ {\psi }_{\mathrm{inside}}\left(x\right)=-c_{\mathrm{left}}(-Q{)}^{-1/4}sin[\frac{1}{\hslash }\underset{x}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}+\frac{\pi }{4}-(\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}+\frac{\pi }{2}\left)\right]\end{array}$

This is identical to the other form of the wave function inside the potential well if the quantity inside curved brackets is an integral multiple of π. The integral between turning points is a nonzero positive quantity. Allowing for a ground state with a quantum number of zero implies choosing integers such that

$\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{-Q}+\frac{\pi }{2}=(n+1)\pi$

which is more nicely written as

$\frac{1}{\hslash }\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{2m(E-V)}=(n+\frac{1}{2})\pi$

This is the WKB quantization condition for a one-dimensional symmetric potential. It is corresponds to Equation (10.5.6) given on page 521 of Bender & Orszag. In addition to this condition one must match coefficients on either side of the potential well:

$-c_{\mathrm{left}}cos\left[\right(n+1\left)\pi \right]=c_{\mathrm{right}}\to c_{\mathrm{left}}=(-1{)}^n{c}_{\mathrm{right}}$

For the Schrödinger equation given at the outset of the presentation, the integral in the quantization condition is

$I=\underset{-x_{\mathrm{turn}}}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{2m(E-\frac{m{\omega }_a^2}{2}|x{|}^a)}=\sqrt{8mE}\underset{0}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{1-\frac{m{\omega }_a^2}{2E}{x}^a}$

For positive exponents, let $X=\frac{m{\omega }_a^2}{2E}{x}^a$ . The resulting integral can be evaluated with a Gauss hypergeometric function:

$\begin{array}{l}I=\sqrt{8mE}(\frac{2E}{m{\omega }_a^2}{)}^{\frac{1}{a}}\frac{1}{a}\underset{0}{\overset{1}{\int }}dX{X}^{\frac{1}{a}-1}\sqrt{1-X}\\ \phantom{I}=\sqrt{8mE}(\frac{2E}{m{\omega }_a^2}{)}^{\frac{1}{a}}{X}^{\frac{1}{a}}{}_{2}F_1[-\frac{1}{2},\frac{1}{a},\frac{1}{a}+1,X]{|}_0^1\\ I=\sqrt{8mE}(\frac{2E}{m{\omega }_a^2}{)}^{\frac{1}{a}}{}_{2}F_1[-\frac{1}{2},\frac{1}{a},\frac{1}{a}+1,1]\end{array}$

The integration is most conveniently executed in more general form to avoid the two-term evaluation returned by Mathematica, then followed by appropriate substitutions. The quantization condition is thus

$\sqrt{8mE}(\frac{2E}{m{\omega }_a^2}{)}^{\frac{1}{a}}{}_{2}F_1[-\frac{1}{2},\frac{1}{a},\frac{1}{a}+1,1]=\hslash (n+\frac{1}{2})\pi$

Evaluating the hypergeometric function using the Gauss summation formula

${}_{2}F_1(a,b;c;1)=\frac{\Gamma \left(c\right)\Gamma (c-a-b)}{\Gamma (c-a)\Gamma (c-b)}$

and solving for energy this becomes

$E=[\hslash \sqrt{\frac{\pi }{2m}}(\frac{m{\omega }_a^2}{2}{)}^{\frac{1}{a}}\frac{\Gamma (\frac{1}{a}+\frac{3}{2})}{\Gamma (\frac{1}{a}+1)}(n+\frac{1}{2}){]}^{\frac{2a}{a+2}},a>0$

For the simple harmonic oscillator a = 2 this expression becomes

$E_{\mathrm{SHO}}=\hslash {\omega }_2(n+\frac{1}{2}),a=2$

which is exactly the same as the analytic eigenvalue solution.

The gamma functions in the energy eigenvalues for arbitrary exponent could be simplified further, but for purposes of visualization it makes no difference. For simplicity of code, let $\hslash =m={\omega }_a=1$ as in the previous presentation. The energy eigenvalue curves for positive exponents are

This graphic is to be compared the corresponding graphic from the previous presentation, for which the eigenvalue data is precalculated and embedded in the page:

Approximations generated by WKB theory are only expected to be accurate for large quantum numbers, but these two graphs look satisfyingly similar even for small quantum numbers. They are compared, one level at a time, in the following interactive graphic that uses data from the previous one for thirty basis states. The WKB approximation is the red line:

As expected the WKB approximation is not accurate for the lowest state, but becomes surprisingly accurate for the other states over most of the domain of the continuous exponent. The only major discrepancy is in the region approaching a = 0 , where the potential becomes constant and there is no quantization. In this limit the WKB approximation to the energy becomes

$\begin{array}{l}\underset{a\to 0}{lim}[\hslash \sqrt{\frac{\pi }{2m}}(\frac{m{\omega }_a^2}{2}{)}^{\frac{1}{a}}\frac{(\frac{1}{a}+\frac{3}{2}{)}^{\frac{1}{a}+1}{e}^{-\frac{1}{a}-\frac{3}{2}}}{(\frac{1}{a}+1{)}^{\frac{1}{a}+\frac{1}{2}}{e}^{-\frac{1}{a}-1}}(n+\frac{1}{2}){]}^{\frac{2a}{a+2}}\\ =\underset{a\to 0}{lim}(\frac{m{\omega }_a^2}{2}\frac{\frac{1}{a}+\frac{3}{2}}{\frac{1}{a}+1}{)}^{\frac{2}{a+2}}[\hslash \sqrt{\frac{\pi }{2m}}\frac{(\frac{1}{a}+\frac{3}{2})}{(\frac{1}{a}+1{)}^{\frac{1}{2}}}{e}^{-\frac{1}{2}}(n+\frac{1}{2}){]}^{\frac{2a}{a+2}}\\ \hfill =\frac{m{\omega }_0^2}{2}\times \underset{a\to 0}{lim}{a}^{-a/(a+2)}=\frac{m{\omega }_0^2}{2}\end{array}$

using an approximation to the gamma function for large argument. Since there is no motion for a constant potential, the total energy is equal to the value of the constant energy. The WKB approximation thus approaches the classical value of the energy for an exponent of zero, which with the chosen constants is one half.

Having an easily computable expression for approximate energy eigenvalues for arbitrary exponent, the wave functions would at first seem to be a simple matter of numerical integration. Since the system is symmetric, a naïve piecewise definition of the wave functions on either side of the origin gives the result

The red line is the potential itself, shifted downward by the energy for the state, so that turning points are where it crosses the x-axis.

Clearly this WKB approximation breaks down at the classical turning points, seen here as truncated spikes on either side of the wave functions. This occurs because the Airy functions applicable at the turning points have not been included. While this could be done by adding an additional part to the piecewise function in the region of the turning point, there is a much more elegant solution: the uniform approximation on page 510 of Bender & Orszag. For wave functions normalized to pass near unity in the ground state, the approximation is

$\psi \left(x\right)=\sqrt{\pi }\left[Q\right(x){]}^{-1/4}(\frac{3}{2}\underset{x_{\mathrm{turn}}}{\overset{x}{\int }}\sqrt{Q\left(x\right)}{)}^{1/6}Ai\left[(\frac{3}{2}\underset{x_{\mathrm{turn}}}{\overset{x}{\int }}\sqrt{Q\left(x\right)}{)}^{2/3}\right]$

As explained in the reference, this single analytic expression covers the regions on either side of the turning point as well as the turning point itself. One problem in a numerical implementation is that complex arguments can get lost because the default range ( −π, π ) is numerically indistinguishable from other complex sheets. One must instead adjust the sign of Q inside the turning points to remove the imaginary unit, as can be seen in the underlying code in the source of this page. It also turns out to be necessary to give the argument of the Airy function an explicit negative sign inside the turning points to ensure a sinusoidal form, with an accompanying reverse in the direction of integration.

The resulting uniform approximation to the wave functions for positive exponents is

Et voila! The truncated spikes at the turning points are gone.

In order to extend the discussion to negative values for the exponent of the potential, the Schrödinger equation given at the outset must be slightly modified. Since the potential for negative exponents does not have a potential well in which energy levels can be quantized, the potential must be inverted with a negative coefficient. The resulting potential well will then have a singularity at the origin, but this does not prevent evaluation of energy levels and wave functions.

The modified Schrödinger equation for negative exponents is

$-\frac{{\hslash }^2}{2m}\frac{d^2\psi }{d{x}^2}-\frac{m{\omega }_a^2}{2}|x{|}^{-\left|a\right|}\psi =E\psi$

with an explicit negative sign in the exponent. Energy eigenvalues are also negative, so that explicit signs will help track which quantities are positive. The integral in the quantization condition is

$I=\underset{0}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{2m(E-V)}=\sqrt{8m(-E)}\underset{0}{\overset{x_{\mathrm{turn}}}{\int }}dx\sqrt{\frac{m{\omega }_a^2}{2(-E)}{x}^{-\left|a\right|}-1}$

where an explicit negative sign has been introduced. Let $\frac{1}{X}=\frac{m{\omega }_a^2}{2(-E)}{x}^{-\left|a\right|}$ and rearrange:

$\begin{array}{l}I=\sqrt{8m(-E)}(\frac{m{\omega }_a^2}{2(-E)}{)}^{\frac{1}{\left|a\right|}}\frac{1}{\left|a\right|}\underset{0}{\overset{1}{\int }}dX{X}^{\frac{1}{\left|a\right|}-1}\sqrt{\frac{1}{X}-1}\\ \phantom{I}=\sqrt{8m(-E)}(\frac{m{\omega }_a^2}{2(-E)}{)}^{\frac{1}{\left|a\right|}}\frac{1}{\left|a\right|}\underset{0}{\overset{1}{\int }}dX{X}^{\frac{1}{\left|a\right|}-\frac{3}{2}}\sqrt{1-X}\\ I=\sqrt{8m}(-E{)}^{\frac{1}{2}-\frac{1}{\left|a\right|}}(\frac{m{\omega }_a^2}{2}{)}^{\frac{1}{\left|a\right|}}\frac{\frac{1}{\left|a\right|}}{\frac{1}{\left|a\right|}-\frac{1}{2}}{}_{2}F_1[-\frac{1}{2},\frac{1}{\left|a\right|}-\frac{1}{2},\frac{1}{\left|a\right|}+\frac{1}{2},1]\end{array}$

Evaluating the hypergeometric function with the Gauss summation formula and inserting in the quantization condition, the resulting energy curves for negative exponents are

$E=-[\hslash \sqrt{\frac{\pi }{2m}}(\frac{2}{m{\omega }_a^2}{)}^{\frac{1}{\left|a\right|}}\frac{\Gamma \left(\frac{1}{\left|a\right|}\right)}{\Gamma (\frac{1}{\left|a\right|}-\frac{1}{2})}(n+\frac{1}{2}){]}^{\frac{2\left|a\right|}{\left|a\right|-2}},-2<a<0$

For the Coulomb-like inverse potential a = −1 this expression is

$E_{\mathrm{inverse}}=-\frac{m^3{\omega }_{-1}^4}{2{\hslash }^2(n+\frac{1}{2}{)}^2}$

The offset of one half to the quantum numbers is presumably an artifact of the WKB approximation, but makes the result curiously proportional to the eigenvalues of a two-dimensional analog of the hydrogen atom.

Again setting $\hslash =m={\omega }_a=1$ , the energy eigenvalue curves for negative exponents greater than minus two are

and the uniform approximation to the wave functions is

The red line is again the potential itself, shifted downward by the energy for the state.

In the limit a → 0 the energy eigenvalue curves are expected to approach the same absolute value with a different sign. This is easily verified for the negative of the energy expression for negative exponents using the same approximation to the gamma function for large argument:

$\begin{array}{l}\underset{a\to 0}{lim}[\hslash \sqrt{\frac{\pi }{2m}}(\frac{2}{m{\omega }_a^2}{)}^{\frac{1}{\left|a\right|}}\frac{(\frac{1}{\left|a\right|}{)}^{\frac{1}{\left|a\right|}-\frac{1}{2}}{e}^{-\frac{1}{\left|a\right|}}}{(\frac{1}{\left|a\right|}-\frac{1}{2}{)}^{\frac{1}{\left|a\right|}-1}{e}^{-\frac{1}{\left|a\right|}+\frac{1}{2}}}(n+\frac{1}{2}){]}^{\frac{2\left|a\right|}{\left|a\right|-2}}\\ =\underset{a\to 0}{lim}(\frac{2}{m{\omega }_a^2}\frac{\frac{1}{\left|a\right|}}{\frac{1}{\left|a\right|}-\frac{1}{2}}{)}^{\frac{2}{\left|a\right|-2}}[\hslash \sqrt{\frac{\pi }{2m}}\frac{(\frac{1}{\left|a\right|}-\frac{1}{2})}{(\frac{1}{\left|a\right|}{)}^{\frac{1}{2}}}{e}^{-\frac{1}{2}}(n+\frac{1}{2}){]}^{\frac{2\left|a\right|}{\left|a\right|-2}}\\ \hfill =\frac{m{\omega }_0^2}{2}\times \underset{a\to 0}{lim}|a{|}^{-|a|/(|a|-2)}=\frac{m{\omega }_0^2}{2}\end{array}$

Given this equal limit, one is tempted to plot the energy eigenvalue curves together with a sign change for negative exponents. The resulting colorful graph

shows which quantum numbers connect across the origin in exponent space. If this plot is magnified the connection appears quite smooth, so that one is then tempted to look for an analytic continuation between the two sets of curves, which unfortunately goes awry. With an identity for the gamma function

$\Gamma (1-x)\Gamma \left(x\right)=\frac{\pi }{sin\pi x}$

rewrite the ratio of gamma functions in the expression for positive energy eigenvalue curves with an explicit negative sign for the exponent:

$\frac{\Gamma (-\frac{1}{\left|a\right|}+\frac{3}{2})}{\Gamma (-\frac{1}{\left|a\right|}+1)}=\frac{\Gamma \left(\frac{1}{\left|a\right|}\right)sin\frac{\pi }{\left|a\right|}}{\Gamma (\frac{1}{\left|a\right|}-\frac{1}{2})sin(\frac{\pi }{\left|a\right|}-\frac{\pi }{2})}=-\frac{\Gamma \left(\frac{1}{\left|a\right|}\right)}{\Gamma (\frac{1}{\left|a\right|}-\frac{1}{2})}tan\frac{\pi }{\left|a\right|}$

The ratio of gamma functions here is that required for the expression for negative exponents. The tangent function is at first unexpected until one realizes that the numerator on the left-hand side is singular whenever the exponent is two divided by an odd integer greater than one. The lack of simple continuation between the two regions in exponent space occurs because these singularities must be removed explicitly.

It is worth remembering that classical orbits for power potentials have the center of force at the center of an ellipse for positive exponents, but at a focus of an ellipse for negative exponents greater than minus two. The lack of simple continuation in the energy eigenvalue expressions represents a quantum mechanical analog of this difference in the two exponent domains.

Uploaded 2014.11.27 — Updated 2017.08.05 analyticphysics.com